Binary effects of sediment on the survival and settlement of larval corals

Here I will explore how (1) DEPOSITED AND (2) SUSPENDED sediment affects the SETTLEMENT and (3) SUSPENDED sediment affects the SURVIVAL of tropical, scleractinean coral LARVAE. The database to be used is comprised of data extracted from studies deemed relevant during the systematic literature review process.

The Dataset

Deposited Sediment

Specifically, we will explore the results of 54 datapoints from 4 studies conducted in 1 ocean on 2 species within 2 genera.

##    n
## 1 54

##    Ref  n
## 1 DS19 38
## 2 DS01  8
## 3 DS06  6
## 4 DS02  2

##                    Ref_name  n
## 1     Ricardo et al. (2017) 38
## 2 Babcock and Davies (1991)  8
## 3           Hodgson (1990a)  6
## 4  Babcock and Smith (2000)  2

##     Ocean  n
## 1 Pacific 54

##                      Gsp  n
## 1     Acropora_millepora 48
## 2 Pocillopora_damicornis  6

##   Updated_Genus  n
## 1      Acropora 48
## 2   Pocillopora  6

Suspended Sediment

For larval settlement, we will explore the results of 30 datapoints from 7 studies conducted in 2 oceans on 4 species within 3 genera.

##    n
## 1 30

##     Ref n
## 1 SS13b 5
## 2 SS13c 5
## 3 SS13d 5
## 4 SS23a 4
## 5 SS24a 4
## 6 SS24b 4
## 7 SS11c 3

##                  Ref_name  n
## 1    Humanes et al. 2017b 15
## 2                 Te 1992  8
## 3 Rushmore 2016 Chapter 2  4
## 4            Gilmour 1999  3

##      Ocean  n
## 1  Pacific 26
## 2 Atlantic  4

##                      Gsp  n
## 1        Acropora_tenuis 15
## 2 Pocillopora_damicornis  8
## 3     Porites_astreoides  4
## 4    Acropora_digitifera  3

##   Updated_Genus  n
## 1      Acropora 18
## 2   Pocillopora  8
## 3       Porites  4

For larval survival, we will explore the results of 63 datapoints from 4 studies conducted in 1 ocean on 5 species within 2 genera.

##    n
## 1 63

##     Ref  n
## 1 SS20d 37
## 2 SS13b  5
## 3 SS13c  5
## 4 SS13d  5
## 5 SS24a  4
## 6 SS24b  4
## 7 SS11b  3

##               Ref_name  n
## 1  Ricardo et al. 2016 37
## 2 Humanes et al. 2017b 15
## 3              Te 1992  8
## 4         Gilmour 1999  3

##     Ocean  n
## 1 Pacific 63

##                      Gsp  n
## 1        Acropora_tenuis 31
## 2     Acropora_millepora 16
## 3 Pocillopora_damicornis  8
## 4      Pocillopora_acuta  5
## 5    Acropora_digitifera  3

##   Updated_Genus  n
## 1      Acropora 50
## 2   Pocillopora 13

Exploratory plots

First let’s explore all data from all species for which there is data about ‘limited settlement’ and ‘reduced larval survival’ as a result of exposure to sediment.

DEFINITIONS:

• ‘Limited settlement’ is defined as a reduction in the number of larval corals settling/attaching to a surface with sediment, as compared to surfaces without sediment (or “control” conditions).
  
• ‘Reduced larval survival’ is defined as a reduction in the number of pre-settlement, larval corals that survive exposure to sediment, as compared to conditions without sediment (or “control” conditions).

Threshold Analyses with Binary Data

Now let’s calculate the thresholds, based on the binary data explored above.

DEFINITIONS:

• ‘LOAEL’, or the ‘lowest observed adverse effect level’ is defined as the lowest dose at which there was an observed adverse effect.
  
• ‘NOAEL’, or the ‘no observed adverse effect level’ is defined as the highest dose at which there was NOT an observed adverse effect.
  
• ‘Adverse effect’ is defined here as any response of a coral individual, colony, or experimental treatment group that may negatively affect the coral’s fitness and/or survival. These adverse effects may include sublethal physiological changes (e.g., significantly reduced growth or photosynthetic rates when compared to coral in ambient/control conditions), bleaching/tissue loss, and mortality. This definition of adverse effect is independent of its magnitude, so while the affect may have potentially reduced fitness, its magnitude may be sufficiently small that the fitness effect is not measureable.

DEPOSITED SEDIMENT

Limited settlement

NOAEL/LOAEL defined by exposure concentration

##   LOAEL
## 1     1

##   NOAEL
## 1     1

SUSPENDED SEDIMENT

Limited settlement

NOAEL/LOAEL defined by exposure concentration

##   LOAEL
## 1 57.79

##   NOAEL
## 1  34.6

Larval mortality

NOAEL/LOAEL defined by exposure concentration

##   LOAEL
## 1    30

##   NOAEL
## 1  29.5

Logistic Regression Model Fitting

Deposited Sediment

##    n
## 1 45

## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular

## Data: settleDS2
## Models:
## modDS4: Binary_limited_settlement ~ Sed_level_stand_mg + (1 | Gsp)
## modDS7: Binary_limited_settlement ~ Sed_level_stand_mg + (1 | Ref)
## modDS10: Binary_limited_settlement ~ Sed_level_stand_mg + (1 | Gsp) + 
## modDS10:     (1 | Ref)
##         npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## modDS4     3 41.580 47.000 -17.790   35.580                         
## modDS7     3 41.515 46.935 -17.758   35.515 0.0651  0     <2e-16 ***
## modDS10    4 43.515 50.742 -17.758   35.515 0.0000  1          1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##        df      AIC
## modDS1  2 39.58049
## modDS7  3 41.51537

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Binary_limited_settlement ~ Sed_level_stand_mg + (1 | Ref)
##    Data: settleDS2
## 
##      AIC      BIC   logLik deviance df.resid 
##     41.5     46.9    -17.8     35.5       42 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7321 -0.7847  0.0016  0.4885  1.2405 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Ref    (Intercept) 1.089    1.044   
## Number of obs: 45, groups:  Ref, 4
## 
## Fixed effects:
##                    Estimate Std. Error z value Pr(>|z|)  
## (Intercept)        -0.06776    0.94288  -0.072   0.9427  
## Sed_level_stand_mg  0.07684    0.04664   1.648   0.0994 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Effect estimates

## [1] 0.8444444

##    
## p    0  1
##   0  9  4
##   1  3 29

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

## Area under the curve: 0.9167

##                   R2m       R2c
## theoretical 0.9847224 0.9885223
## delta       0.9765060 0.9802742

## $Ref

##  Groups Name        Std.Dev.
##  Ref    (Intercept) 1.0436

## [1] 2.839528

##                          Est        LL       UL
## (Intercept)        0.9541154 0.2650232 3.434930
## Sed_level_stand_mg 1.0547038 0.9899510 1.123692

There is suggestive, but non-significant, evidence that for every doubling of exposure concentration of deposited sediment, the odds of reduced settlement rate increase by 1.05 times (95% CI 0.99, 1.12; GLMM z = 1.648, p = 0.09) after accounting for variability among studies.

Prediction figure

#Overlaying glmm results on figure, but prediction line is jagged...
pred_modDS7 <- predict(modDS7, newdata=settleDS2, type="response")
settleDS3 <- cbind(settleDS2, pred_modDS7)

settleDS_plot <- ggplot(data = settleDS3) +
    geom_point(mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_limited_settlement,
      color = Ref)) +
    labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"), 
         y = "Reduced settlement rate due to sediment exposure",
         color = "Study") +
    scale_x_log10(limits=c(0.1,1100),breaks=c(0.01,0.1,1,10,100,1000),
                  label=c("0.01","0.1","1","10","100","1000")) +
    geom_line(aes(x = Sed_level_stand_mg, y = pred_modDS7), inherit.aes=FALSE) +
    theme_bw()

settleDS_plot

That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/

jvalues <- with(settleDS2, seq(from = min(Sed_level_stand_mg), to = max(Sed_level_stand_mg), length.out = 100))

# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
    settleDS2$Sed_level_stand_mg <- j
    predict(modDS7, newdata = settleDS2, type = "response")
})

# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
    c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))

# add in Sed_level_stand_mg values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))

# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability", 
                       "LowerQuantile", "UpperQuantile", "Sed_conc")
#head(plotdat)

# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_conc)) + 
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) +
  ylim(c(0, 1))

ggplot(plotdat, aes(x = Sed_conc)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1))

#Overlaying average marginal predicted probabilities on figure
#by Ref
settleDS_plot3 <- ggplot() +
    geom_point(data = settleDS2, mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_limited_settlement,
      color = Ref)) +
    labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"), 
         y = "Predicted probability of reduced settlement\nrate due to sediment exposure",
         color = "Binary Data\nby Study",
         linetype = "Predicted\nProbability") +
    scale_x_log10(limits=c(0.1,1100), breaks=c(0.01,0.1,1,10,100,1000,loaelDS),
                  label=c("0.01","0.1","1","10","100","1000",round(loaelDS, digits=1))) +
    geom_ribbon(data = plotdat, aes(x = Sed_conc, y = PredictedProbability, 
                                    ymin = LowerQuantile, ymax = UpperQuantile), 
                alpha = 0.15) +
    geom_line(data = plotdat, aes(x = Sed_conc, y = PredictedProbability, 
                                  linetype = "twodash")) +
    geom_line(data = plotdat, aes(x = Sed_conc, y = MedianProbability, 
                                  linetype = "solid")) +
    geom_vline(xintercept=loaelDS, linetype="dashed", color = "red") +
    theme_bw() +
    scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))

settleDS_plot3

Suspended Sediment

##    n
## 1 20

## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular

## Data: settleSS2
## Models:
## modSS4: Binary_limited_settlement ~ Sed_level_stand_mg + (1 | Gsp)
## modSS7: Binary_limited_settlement ~ Sed_level_stand_mg + (1 | Ref)
## modSS10: Binary_limited_settlement ~ Sed_level_stand_mg + (1 | Gsp) + 
## modSS10:     (1 | Ref)
## modSS13: Binary_limited_settlement ~ Sed_level_stand_mg + (1 | Ref_name/Ref)
##         npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## modSS4     3 20.127 23.114 -7.0636   14.127                         
## modSS7     3 19.416 22.403 -6.7079   13.416 0.7115  0     <2e-16 ***
## modSS10    4 22.127 26.110 -7.0636   14.127 0.0000  1          1    
## modSS13    4 22.127 26.110 -7.0636   14.127 0.0000  0          1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##        df      AIC
## modSS1  2 20.75610
## modSS7  3 19.41575

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Binary_limited_settlement ~ Sed_level_stand_mg + (1 | Ref)
##    Data: settleSS2
## 
##      AIC      BIC   logLik deviance df.resid 
##     19.4     22.4     -6.7     13.4       17 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -0.55090 -0.01777 -0.01380 -0.01347  1.72158 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Ref    (Intercept) 120.9    10.99   
## Number of obs: 20, groups:  Ref, 6
## 
## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)
## (Intercept)        -8.522305   5.925743  -1.438    0.150
## Sed_level_stand_mg  0.001455   0.005675   0.256    0.798
## convergence code: 0
## Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Binary_limited_settlement ~ log10_conc + (1 | Ref)
##    Data: settleSS2
## 
##      AIC      BIC   logLik deviance df.resid 
##     14.3     17.3     -4.1      8.3       17 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -0.01213  0.00000  0.00000  0.00000  0.05431 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Ref    (Intercept) 18037    134.3   
## Number of obs: 20, groups:  Ref, 6
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   -95.32      46.01  -2.072   0.0383 *
## log10_conc     27.27      14.22   1.918   0.0551 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Effect estimates

## [1] 1

##    
## p    0  1
##   0 17  0
##   1  0  3

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

## Area under the curve: 1

##                      R2m          R2c
## theoretical 1.781484e-02 9.998209e-01
## delta       7.265477e-14 4.077598e-12

## $Ref

##  Groups Name        Std.Dev.
##  Ref    (Intercept) 134.3

## [1] 2.118497e+58

##                      Est           LL           UL
## (Intercept) 4.001764e-42 2.730730e-81 5.864407e-03
## log10_conc  6.978395e+11 5.493766e-01 8.864228e+23

There is suggestive, though non-significant, evidence that for every 10-fold increase in exposure concentration of suspended sediment, the odds of reduced settlement rate increase (GLMM z = 1.918, p = 0.055), while accounting for variability among studies.

Prediction figures

#Overlaying glmm results on figure, but prediction line is jagged...
pred_modSS7_log <- predict(modSS7_log, newdata=settleSS2, type="response")
settleSS3 <- cbind(settleSS2, pred_modSS7_log)

settleSS_plot <- ggplot(data = settleSS3) +
    geom_point(mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_limited_settlement,
      color = Ref)) +
    labs(x = "Sediment exposure concentration (mg/L)", 
         y = "Reduced settlement rate\ndue to sediment exposure",
         color = "Study") +
    scale_x_log10(limits=c(1,1000),breaks=c(0.01,0.1,1,10,100,1000),
                  label=c("0.01","0.1","1","10","100","1000")) +
    geom_line(aes(x = Sed_level_stand_mg, y = pred_modSS7_log), inherit.aes=FALSE) +
    theme_bw()

settleSS_plot

That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/

By sediment exposure concentration

jvalues <- with(settleSS2, seq(from = min(log10_conc), to = max(log10_conc), length.out = 100))

# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
    settleSS2$log10_conc <- j
    predict(modSS7_log, newdata = settleSS2, type = "response")
})

# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
    c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))

# add in log10_conc values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
plotdat <- plotdat %>% mutate(Sed_level_linear=10^jvalues)

# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability", 
                       "LowerQuantile", "UpperQuantile", "log10_mg_L", "mg_L")
#head(plotdat)

# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = log10_mg_L)) + 
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) +
  ylim(c(0, 1))

ggplot(plotdat, aes(x = log10_mg_L)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1))

ggplot(plotdat, aes(x = mg_L, y = PredictedProbability)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1)) +
  scale_x_log10()

#Overlaying average marginal predicted probabilities on figure
#by Ref
settleSS_plot2 <- ggplot() +
    geom_point(data = settleSS2, mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_limited_settlement,
      color = Ref)) +
    labs(x = "Sediment exposure concentration (mg/L)", 
         y = "Predicted probability of reduced settlement\nrate due to sediment exposure",
         color = "Binary Data\nby Study",
         linetype = "Predicted\nProbability") +
    scale_x_log10(limits=c(1,max(settleSS2$Sed_level_stand_mg)), 
                  breaks=c(0.01,0.1,1,10,100,1000,loaelSS1),
                  label=c("0.01","0.1","1","10","100","1000",round(loaelSS1,digits = 1))) +
    geom_ribbon(data = plotdat, aes(x = mg_L, y = PredictedProbability, 
                                    ymin = LowerQuantile, ymax = UpperQuantile), 
                alpha = 0.15) +
    geom_line(data = plotdat, aes(x = mg_L, y = PredictedProbability, 
                                  linetype = "twodash")) +
    geom_line(data = plotdat, aes(x = mg_L, y = MedianProbability, 
                                  linetype = "solid")) +
    geom_vline(xintercept=loaelSS1, linetype="dashed", color = "red") +
    theme_bw() +
    scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))

settleSS_plot2

Larval Survival

Suspended Sediment

##    n
## 1 52

##          df      AIC
## modSSb1   2 44.75628
## modSSb1b  3 40.72043

## Data: survivalSS2
## Models:
## modSSb4: Binary_larval_survival ~ Sed_level_stand_mg + (1 | Gsp)
## modSSb7: Binary_larval_survival ~ Sed_level_stand_mg + (1 | Ref)
## modSSb10: Binary_larval_survival ~ Sed_level_stand_mg + (1 | Gsp) + (1 | 
## modSSb10:     Ref)
## modSSb13: Binary_larval_survival ~ Sed_level_stand_mg + (1 | Ref_name/Ref)
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## modSSb4     3 44.117 49.971 -19.059   38.117                         
## modSSb7     3 43.904 49.758 -18.952   37.904 0.2132  0     <2e-16 ***
## modSSb10    4 44.063 51.868 -18.031   36.063 1.8415  1     0.1748    
## modSSb13    4 45.695 53.500 -18.848   37.695 0.0000  0     1.0000    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##         df      AIC
## modSSb1  2 44.75628
## modSSb4  3 44.11739

## 
## Call:
## glm(formula = Binary_larval_survival ~ Sed_level_stand_mg, family = binomial(link = "logit"), 
##     data = survivalSS2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.7033  -0.5272  -0.5016  -0.4975   2.0679  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        -2.0356448  0.5215749  -3.903 9.51e-05 ***
## Sed_level_stand_mg  0.0007646  0.0012874   0.594    0.553    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 41.087  on 51  degrees of freedom
## Residual deviance: 40.756  on 50  degrees of freedom
## AIC: 44.756
## 
## Number of Fisher Scoring iterations: 4

Effect estimates

## [1] 0.8653846

##    
## p    0  1
##   0 45  7

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

## Area under the curve: 0.6429

##                     R2m         R2c
## theoretical 0.014274072 0.014274072
## delta       0.005519128 0.005519128

There is no significant relationship between exposure concentration of suspended sediment and the odds of larval mortality (GLMM z = 0.762, p = 0.446).

Prediction figures

#Overlaying glmm results on figure, but prediction line is jagged...
pred_modSSb1 <- predict(modSSb1, newdata=survivalSS2, type="response")
survivalSS3 <- cbind(survivalSS2, pred_modSSb1)

survivalSS_plot <- ggplot(data = survivalSS3) +
    geom_point(mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_larval_survival,
      color = Ref)) +
    labs(x = "Sediment exposure concentration (mg/L)", 
         y = "Larval mortality due to sediment exposure",
         color = "Study") +
    scale_x_log10(limits=c(1,1000),breaks=c(0.01,0.1,1,10,100,1000),
                  label=c("0.01","0.1","1","10","100","1000")) +
    geom_line(aes(x = Sed_level_stand_mg, y = pred_modSSb1), inherit.aes=FALSE) +
    theme_bw()

survivalSS_plot

That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/

By sediment exposure concentration

jvalues <- with(survivalSS2, seq(from = min(Sed_level_stand_mg), to = max(Sed_level_stand_mg), length.out = 100))

# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
    survivalSS2$Sed_level_stand_mg <- j
    predict(modSSb1, newdata = survivalSS2, type = "response")
})

# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
    c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))

# add in Sed_level_stand_mg values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))

# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability", 
                       "LowerQuantile", "UpperQuantile", "Sed_conc")
#head(plotdat)

# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_conc)) + 
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) +
  ylim(c(0, 1))

ggplot(plotdat, aes(x = Sed_conc)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1))

#Overlaying average marginal predicted probabilities on figure
#by Ref
survivalSS_plot2 <- ggplot() +
    geom_point(data = survivalSS2, mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_larval_survival,
      color = Ref)) +
    labs(x = "Sediment exposure concentration (mg/L)", 
         y = "Predicted probability of larval\nmortality due to sediment exposure",
         color = "Binary Data\nby Study",
         linetype = "Predicted\nProbability") +
    scale_x_log10(limits=c(1,max(survivalSS2$Sed_level_stand_mg)), 
                  breaks=c(0.01,0.1,1,10,100,1000,loaelSS2),
                  label=c("0.01","0.1","1","10","100","1000",round(loaelSS2,digits = 1))) +
    geom_ribbon(data = plotdat, aes(x = Sed_conc, y = PredictedProbability, 
                                    ymin = LowerQuantile, ymax = UpperQuantile), 
                alpha = 0.15) +
    geom_line(data = plotdat, aes(x = Sed_conc, y = PredictedProbability, 
                                  linetype = "twodash")) +
    geom_line(data = plotdat, aes(x = Sed_conc, y = MedianProbability, 
                                  linetype = "solid")) +
    geom_vline(xintercept=loaelSS2, linetype="dashed", color = "red") +
    theme_bw() +
    scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))

survivalSS_plot2